\chapter{Analisi teorica}

In questa sezione ci occuperemo di descrivere gli aspetti teorici relativi ai problemi differenziali che intendiamo risolvere per mezzo della nostra libreria.
Porremo particolare attenzione all'analisi dei metodi spettrali, in quanto rappresentano il vero fulcro attorno al quale ruota tutta la nostra trattazione.
Infatti, in questo progetto, abbiamo scelto di risolvere problemi differenziali evolutivi per mezzo di metodi spettrali che, a differenza dei classici metodi agli elementi finiti, hanno il vantaggio di essere di grado più elevato.
In particolare, la discretizzazione spaziale di un problema ellittico attraverso il metodo degli elementi finiti (FEM) vincola l'ordine di convergenza della soluzione numerica al grado dei polinomi usati.
Invece la discretizzazione spaziale attraverso metodi spettrali (SM) e metodi agli elementi spettrali (SEM) garantisce una velocità di convergenza ottimale, limitata solo dalla regolarità della soluzione del problema.
In prima approssimazione, l'idea di base di questi metodi è quella di utilizzare polinomi \emph{globali} sul dominio computazionale, invece di polinomi \emph{a tratti} sui singoli elementi.


\section{I metodi spettrali}

Ricordiamo innanzitutto i principali strumenti matematici necessari per un inquadramento generale dei metodi spettrali applicati alla risoluzione di problemi ellittici.
In effetti il nostro obiettivo finale sarà quello di risolvere problemi evolutivi ma, di fatto, questo scopo sarà raggiunto attraverso una discretizzazione dell'intervallo temporale e, ad ogni istante, quello che dovremo risolvere sarà comunque un problema ellittico.
Per questo motivo, il seguito del discorso può essere affrontato senza alcun pericolo di sviare dall'obiettivo principale di questo lavoro.\\
Sia $\varOmega$ un insieme aperto in $\R^2$ e $\vec{x}$ un punto di $\varOmega$. Lo spazio $L^2(\varOmega)$ delle funzioni quadrato-sommabili è definito da
\[
	L^2(\varOmega) = \left\{ f: \varOmega \rightarrow \R^2 : \int_{\varOmega} \vert f(\vec{x}) \vert^2 \,\D\varOmega < \infty \right\},
\]
in modo tale che la norma a lui associata sia definita come
\[
	\Vert f \Vert_{L^2(\varOmega)} = \left(\int_{\varOmega} \vert f(\vec{x}) \vert^2 \,\D\varOmega\right)^{1/2}.
\]
Inoltre, definiamo anche lo spazio delle funzioni superiormente limitate come
\[
	L^{\infty}(\varOmega) = \{ f: \varOmega \rightarrow \R\quad \mathrm{limitate\ a\ meno\ di\ insiemi\ di\ misura\ nulla}\},
\]
dove, in questo caso, la norma associata è descritta da
\[
	\Vert f \Vert_{L^{\infty}(\varOmega)} = \underset{\varOmega}{\esssup} \vert f(\vec{x}) \vert.
\]
Per ultimo, vogliamo definire anche lo spazio
\[
	H^1(\varOmega) = \{ f \in L^2(\varOmega) : \nabla f \in L^2(\varOmega)\},
\]
in quanto sarà di fondamentale importanza nel seguito della trattazione, poiché vedremo che sarà lo spazio ideale in cui ambientare i problemi differenziali che ci interesseranno.
In questo caso la sua norma è definita come
\[
	\Vert f \Vert_{H^1(\varOmega)} = \left(\Vert f \Vert^2_{L^2(\varOmega)} + \Vert \nabla f \Vert^2_{L^2(\varOmega)}\right)^{1/2}.
\]

A questo punto, consideriamo un generico problema ellittico bidimensionale di diffusione-reazione con condizioni al contorno miste:
\begin{equation}\label{eq:ellittico}
	\left\{
  	\begin{array}{rll}
  		Lu 					&:= -\Div(\mu \nabla u) + \sigma u = f		&\quad\mathrm{in}\,\varOmega\\[2mm]
  		u 					&= g 								&\quad\mathrm{su}\,\varGamma_D\\[2mm]
  		\nabla u \cdot\vec{n} 	&= h 								&\quad\mathrm{su}\,\varGamma_N
	\end{array}
	\right.
\end{equation}
dove il dominio è $\varOmega \subset \R^2$, con frontiera $\partial\varOmega = \varGamma_D \cup \varGamma_N$, e $\varGamma_D \cap \varGamma_N = \emptyset$.
È bene ricordare che con $\vec{n}$ intendiamo il versore normale uscente dal dominio.\\
La formulazione debole del problema (\ref{eq:ellittico}) è:
\begin{equation}\label{eq:debole}
	\mathrm{Trovare}\ u \in V \quad : \quad a(u,v) = F(v), \qquad \forall v \in V,
\end{equation}
dove $V$ è lo spazio di Hilbert $H^1_{\varGamma_D}(\varOmega)$, definito come
\[
	H^1_{\varGamma_D}(\varOmega) = \left\{ u \in H^1(\varOmega) : u\Big\vert_{\varGamma_D} = 0 \right\}.
\]
Secondo questo formalismo, possiamo dunque definire le forme variazionali:
\begin{eqnarray*}
	a(u,v)	&=	&\int_{\varOmega} \mu \nabla u \nabla v + \int_{\varOmega} \sigma uv\\[2mm]
	F(v)		&=	&\int_{\varOmega}fv + \int_{\varGamma_N} hv.
\end{eqnarray*}
Naturalmente questa formulazione è valida solo nel caso in cui $\varGamma_D = \emptyset$ mentre, nel caso non omogeneo, sarà necessario provvedere a un rilevamento del dato Dirichlet.
Non è obiettivo di questo lavoro scendere in tutti i dettagli riguardanti la formulazione debole dei problemi ellittici ma, in ogni caso, per tutti i dettagli rimandiamo a~\cite{QUARTERONI}.


\subsection{I metodi SM e SEM}

Passiamo ora ad analizzare nel dettaglio il meccanismo di funzionamento vero e proprio dei metodi spettrali.\\
Introduciamo, per ogni numero intero $N\in\N$, lo spazio $\Q_N$ dei polinomi a coefficienti reali, di grado minore o uguale a $N$ rispetto a \emph{ciascuna} delle variabili.
In particolare, nel caso bidimensionale avremo:
\[
	\Q_N (\varOmega) = \left\{ v(\vec{x}) = \sum_{k,m=0}^N a_{km}\,x_1^k x_2^m, \qquad\mathrm{con}\ a_{km} \in \R \right\}.
\]
Nel caso del metodo di Galerkin spettrale SM, lo spazio $V$ della formulazione variazionale viene approssimato con uno spazio a dimensione finita $V_N \subset \Q_N$.
Dato che $V$ è un sottospazio di $H^1_{\varGamma_D}(\varOmega)$, $V_N$ viene scelto come l'insieme dei polinomi di $\Q_N$ che si annullano sulla parte di frontiera $\varGamma_D$, caratterizzata da una condizione di tipo Dirichlet.
Dunque, analogamente a quanto definito nella precedente sezione,
\[
	V_N = \left\{ v_N \in \Q_N \quad : \quad v_N \Big\vert_{\varGamma_D} = 0 \right\}.
\]
L'approssimazione SM del problema (\ref{eq:ellittico}) risulta quindi essere:
\begin{equation}\label{eq:SM}
	\mathrm{Trovare}\ u_N \in V_N \quad : \quad a(u_N,v_N) = F(v_N), \qquad \forall v_N \in V_N.
\end{equation}
Ora restringiamoci a un dominio $\varOmega$ costituito dall'unione di quadrilateri $\varOmega_k$, ciascuno dei quali riconducibile al quadrato di riferimento $\hat{\varOmega} = [-1,1]^2$ mediante una trasformazione invertibile $\boldsymbol{\varphi}_k:\hat{\varOmega} \rightarrow \varOmega_k$.\\
Incidentalmente osserviamo che la formulazione SM non può essere usata nel caso in cui il problema richieda che la soluzione si annulli solo su una porzione di lato.
Infatti, in tal caso, un polinomio che rispettasse una simile caratteristica sarebbe necessariamente nullo su tutto il lato.\\
Per ovviare a questa problematica, e dunque utilizzare un metodo decisamente più maneggevole, si ricorre al metodo agli elementi spettrali (SEM).
In questo caso, l'ambiente discreto che approssima lo spazio infinito dimensionale $V$ è
\[
	V_N^C = \left\{ v_N \in C^0(\overline{\varOmega}) \quad : \quad v_N \Big\vert_{\varOmega_k}\!\! \circ\:\boldsymbol{\varphi}_k \in \Q_N (\hat{\varOmega}) \right\}.
\]
Si può pertanto osservare che, una volta scelto questo ambiente, l'approssimazione SEM del problema (\ref{eq:ellittico}) diventa:
\begin{equation}\label{eq:SEM}
	\mathrm{Trovare}\ u_N \in V_N^C \quad : \quad a_C(u_N,v_N) = F_C(v_N), \qquad \forall v_N \in V_N^C,
\end{equation}
dove
\[
	a_C(u_N,v_N) = \sum_k a_{\varOmega_k}(u_N,v_N), \qquad F_C(v_N) = \sum_k F_{\varOmega_k}(v_N),
\]
con $a_{\varOmega_k}(\cdot\,,\cdot)$ e $F_{\varOmega_k}(\cdot)$ restrizioni di $a(\cdot\,,\cdot)$ e di $F(\cdot)$ al sottodominio $\varOmega_k$.\\
Di fatto possiamo constatare come il metodo agli elementi spettrali assomigli molto, almeno nell'idea, a quello più classico agli elementi finiti.
In questo caso il vantaggio, come già discusso in precedenza, riguarda la possibilità di usare dei polinomi di grado più elevato, che consentano una maggiore accuratezza nel calcolo della soluzione, oltre a un più alto ordine di convergenza.


\subsection{Analisi dei metodi SM e SEM}

Per quanto riguarda l'analisi di esistenza e unicità, oltre che di stabilità e convergenza dei metodi SM e SEM, si sfrutta il lemma di Lax-Milgram, comunque applicabile ai metodi di tipo Galerkin.
\begin{lemma}[Lax-Milgram]
	Siano $V$ uno spazio di Hilbert, $a(u,v)$ una forma bilineare su $V$ e $F$ un funzionale lineare e continuo appartenente allo spazio $V'$, duale di $V$. Se la forma bilineare a è:
	\begin{itemize}
		\item continua:
			\[
				\exists M < \infty \quad : \quad \vert a(u,v)\vert \leq M \Vert u \Vert \, \Vert v \Vert , \qquad \forall u,v \in V;
			\]
		\item coerciva:
			\[
				\exists \alpha > 0 \quad : \quad a(v,v) \geq \alpha \Vert v \Vert^2, \qquad \forall v \in V,
			\]
	\end{itemize}
	allora esiste un'unica soluzione $u \in V$. Inoltre vale la seguente stima di stabilità:
	\begin{equation}\label{eq:stab}
		\Vert u \Vert \leq \frac{1}{\alpha} \Vert F \Vert_{V'}
	\end{equation}
\end{lemma}
La dimostrazione di questo lemma si basa su diversi risultati di analisi funzionale, tra cui il teorema di rappresentazione di Riesz.
In ogni caso, poiché non intendiamo scendere nei particolari, per tutti i dettagli tecnici riguardanti il lemma e la sua dimostrazione rimandiamo a~\cite{SALSA}.\\
Venendo invece ad affrontare il problema della consistenza, poiché abbiamo a che fare con un metodo di Galerkin, seppur spettrale, siamo comunque in grado di applicare il lemma di Ceà.
\begin{lemma}[Ceà]
	Sotto le ipotesi del lemma di Lax-Milgram, supponiamo che $u_N$ sia soluzione del problema approssimato e che $u$ sia soluzione di quello originale.
	Allora il metodo di Galerkin è fortemente consistente, ovvero:
	\[
		a(u-u_N,v_N) = 0, \qquad \forall v_N \in V_N.
	\]
	Inoltre, vale la seguente stima:
	\[
		\Vert u-u_N \Vert \leq \frac{M}{\alpha} \inf_{v \in V_N} \Vert u-v \Vert,
	\]
	dove $M$ e $\alpha$ sono, rispettivamente, le costanti di continuità e di coercività che figurano nell'enunciato del lemma di Lax-Milgram.
\end{lemma}
Come dicevamo, dato che i metodi SM e SEM rappresentano un caso particolare del metodo di Galerkin, dall'applicazione diretta dei lemmi precedenti si hanno i risultati di esistenza e unicità della soluzione, oltre a quelli di stabilità e convergenza del problema.\\
Si possono inoltre dimostrare le seguenti \emph{stime a priori} per l'errore che, come abbiamo fatto fin'ora, mostreremo senza scendere nei dettagli.
Per un approfondimento ulteriore si veda~\cite{SPETTRALI}.
\begin{theorem}
	Sia $u \in V$ la soluzione esatta e $u_N$ la corrispondente soluzione approssimata ottenuta con il metodo SM. Si supponga inoltre che $u \in H^{s+1}(\varOmega)$, per qualche $s \geq 0$. Allora vale la seguente stima:
	\[
		\Vert u-u_N \Vert_{H^1(\varOmega)} \leq C_s N^{-s} \Vert u \Vert_{H^{s+1}(\varOmega)},
	\]
	dove $N$ è il grado dei polinomi approssimanti e $C_s$ una costante che non dipende da $N$, ma può dipendere da $s$.\\
	Invece, se $u_N$ è la soluzione ottenuta con il metodo SEM, in tal caso avremo:
	\[
		\Vert u-u_N \Vert_{H^1(\varOmega)} \leq C_s H^{\min (N,s)} N^{-s} \Vert u \Vert_{H^{s+1},(\varOmega)}
	\]
	con $H$ lunghezza massima dei lati dei singoli elementi $\varOmega_k$.
\end{theorem}
Si può notare come, a differenza di un classico metodo agli elementi finiti, in questo caso la velocità di convergenza non è limitata dalla sola scelta del grado dei polinomi utilizzati nella discretizzazione, bensì anche dalla regolarità della soluzione.
Questo aspetto può trasformarsi in un incredibile vantaggio poiché, nel caso in cui la soluzione esatta sia analitica, e cioè valga
\[
	\exists \gamma >0\quad : \quad \Vert u-u_N \Vert_{H^1(\varOmega)} \leq C \exp (-\gamma N),
\]
allora l'ordine di convergenza diventa addirittura esponenziale.
È questo il principale vantaggio che, utilizzando i metodi spettrali, siamo in grado di sfruttare.\\
Tuttavia tali metodi possono essere applicati soltanto a geometrie semplici, come quadrilateri o unioni di quadrilateri, mappabili con una trasformazione invertibile in un quadrato.
Inoltre, nel caso spettrale, la matrice di stiffness è molto meno sparsa rispetto a quella associata al metodo degli elementi finiti, a causa della struttura globale dei polinomi utilizzati.
Questo si traduce in un notevole sforzo computazionale necessario al calcolo della matrice di stiffness e del termine noto, oltre alla risoluzione algebrica del sistema di equazioni associato al metodo.
Per ovviare a questo problema, si procede utilizzando formule di integrazione numerica di tipo gaussiano, che danno origine ai metodi G-NI e SEM-NI.


\section{Approssimazione polinomiale}

In questa sezione ci occuperemo di tutta la teoria riguardante la discretizzazione numerica degli integrali che compaiono nelle formulazioni deboli dei problemi variazionali.
Come accennavamo in precedenza, tali integrazioni numeriche vengono effettuate con metodi altamente efficienti, i quali consentono di ottenere formule di quadratura con grado di esattezza massimale.
Il modo naturale di raggiungere questo obiettivo è quello di effettuare una scelta ben ponderata dei nodi di integrazione, in modo da soddisfare alcuni principi che fra poco descriveremo.
Poiché una delle strade che conduce a ottenere formule di quadratura con grado di precisione massimale è quella dei \emph{polinomi di Legendre}, è da qui che cominciamo la nostra descrizione.

\subsection{Polinomi di Legendre}

Per cominciare, ci restringiamo all'intervallo monodimensionale $(-1,1)$.
Definiamo innanzitutto il prodotto scalare in $L^2(-1,1)$:
\begin{equation}\label{eq:prodscal}
	(f,g) = \int_{-1}^1 f(x)g(x)\,\D x.
\end{equation}
Come accennavamo in precedenza, i primi strumenti da introdurre per la costruzione di opportune formule di integrazione numerica sono i polinomi di Legendre
\[
	L_k(x),\quad\mathrm{per}\ k = 0,1,2,\dots,
\]
definiti come le autofunzioni del problema di Sturm-Liouville
\[
	\dfrac{\D}{\D x}\left((1-x^2)\dfrac{\D}{\D x}L_k(x)\right) + k(k+1)L_k(x)=0.
\]
Ora, se $L_k(x)$ è normalizzato in modo che $L_k(1) = 1$, allora si può definire
\[
	L_k(x) = \frac{1}{2^k} \sum_{j=0}^{[k/2]} (-1)^j \binom{k}{j} \binom{2k-2j}{k} x^{k-2j},\qquad\forall k=0,1,2,\ldots
\]
dove $[k/2]$ è la parte intera di $k/2$.
In alternativa alla precedente formulazione, è anche possibile scrivere una relazione ricorsiva, in modo che risulti:
\begin{equation}
	\left\{
	\begin{array}{rll}
		L_0(x)		&=	&1\\[4.5mm]
		L_1(x)		&=	&x\\[3mm]
		L_{k+1}(x)		&=	&\dfrac{2k+1}{k+1}xL_k(x)- \dfrac{k}{k+1}L_{k-1}(x),\qquad k=1,2,\ldots
	\end{array}
	\right.
\end{equation}
Tra i vari polinomi intercorre una relazione di ortogonalità, ovverosia
\begin{equation}
	(L_k,L_m) = \left\{
	\begin{array}{cl}
		0 						&\qquad\mathrm{se}\ m \neq k\\[3mm]
		\left(k+\dfrac{1}{2}\right)^{-1}	&\qquad\mathrm{se}\ m=k.\\
	\end{array}
	\right.
\end{equation}
Grazie a questa proprietà, i polinomi di Legendre sono linearmente indipendenti e formano una base per lo spazio $L^2(-1,1)$.\\
A questo punto risulta immediato definire lo sviluppo in \emph{serie di Legendre} di una qualsiasi funzione $f \in L^2(-1,1)$:
\begin{equation}
	f(x) = \sum_{k=0}^{\infty} \hat{f}_k L_k(x).
\end{equation}
Ora, sfruttando il seguente risultato:
\begin{eqnarray*}
	(f,L_k)	&=	& \int_{-1}^1 f(x)L_k(x)\,\D x\\[2mm]
			&=	& \int_{-1}^1 \left( \sum_{i=0}^{\infty} \hat{f}_i L_i(x)\right)L_k(x)\,\D x\\[2mm]
			&=	& \sum_{i=0}^{\infty} \left( \int_{-1}^1 L_i(x)L_k(x)\,\D x \right) \hat{f}_i\\[3mm]
			&=	& \hat{f}_k \Vert L_k \Vert_{L^2(-1,1)}^2
\end{eqnarray*}
è possibile scrivere l'espressione dei \emph{coefficienti di Legendre}:
\begin{equation}
	\hat{f}_k = \dfrac{(f,L_k)}{ \Vert L_k \Vert_{L^2(-1,1)}^2} = \left(k+\dfrac{1}{2}\right)\int_{-1}^1 f(x)L_k(x)\,\D x.
\end{equation}
Dalla conoscenza di tali coefficienti risulta quindi immediato rappresentare una qualsiasi funzione di $L^2(-1,1)$ secondo l'opportuna serie di Legendre.
Consideriamo ora la troncata $N$-esima della serie di Legendre di $f$, cioè
\begin{equation}
	f_N(x) = \sum_{k=0}^N \hat{f_k}L_k(x).
\end{equation}
È proprio questa rappresentazione che ci interesserà nel seguito della trattazione e, a tal proposito, è possibile dimostrare alcuni risultati di convergenza.
\begin{theorem}
	Per ogni $f \in L^2(-1,1)$ la sua serie di Legendre converge a $f$ in norma $L^2(-1,1)$:
	\[
		\lim_{N \rightarrow \infty} \Vert f-f_N \Vert_{L^2(-1,1)} = 0.
	\]
	Inoltre, se $f \in H^s(-1,1)$ per qualche $s \geq 1$, allora esiste una costante $C_s>0$, indipendente da $N$, tale che
	\[
		\Vert f-f_N \Vert_{L^2(-1,1)} \leq C_s \left( \frac{1}{N} \right)^s \left\Vert \dfrac{\D^s}{\D x^s}f \right\Vert_{L^2(-1,1)}
	\]
\end{theorem}
Anche in questo caso riportiamo il risultato, fondamentale per garantire la sensatezza dell'approccio utilizzato, evitando di scendere nei dettagli della dimostrazione per non divagare dal discorso principale.
Passiamo invece a descrivere nel dettaglio il funzionamento delle formule di quadratura che saranno utilizzate nella libreria.


\subsection{Formule di quadratura}

Sia $N$ il numero di nodi di quadratura della formula che vogliamo prendere in considerazione.
In realtà la nostra scelta verterà su formule di quadratura interpolatorie del tipo
\begin{equation}
	\int_{-1}^1 f(x)\,\D x = \sum_{i=0}^N \alpha_i f(x_i),
\end{equation}
dunque formule di grado massimo ottenibile usando $N$ nodi $\{ x_i\}$.\\
Introduciamo i nodi di quadratura di Gauss-Legendre-Lobatto (GLL), consistenti negli estremi dell'intervallo $(-1,1)$ e negli zeri della derivata prima del polinomio di Legendre $L_N$.\\
In questo modo i polinomi caratteristici di Lagrange corrispondenti a tali nodi
\[
	\psi_i \in \Q_N \quad : \quad \psi_i(x_j) = \delta_{ij}, \qquad 0 \leq i,j \leq N
\]
hanno forma analitica
\begin{equation}
	\psi_i (x) = \dfrac{-1}{N(N+1)} \dfrac{(1-x^2)}{(x-x_i)L_N(x_i)}\dfrac{\D}{\D x}L_N(x), \qquad i=0,\dots,N.
\end{equation}
Pertanto definiamo, per ogni funzione $f \in C^0([-1,1])$, il suo polinomio d'interpolazione $\varPi_N^{GLL} \in \Q_N$ nei nodi GLL:
\begin{equation}
	\varPi_N^{GLL} f(x) = \sum_{i=0}^N f(x_i) \psi_i(x).
\end{equation}
In base a come è stato definito, per questo polinomio vale la relazione:
\[
	\varPi_N^{GLL}f(x_i) = f(x_i), \qquad 0 \leq i \leq N,
\]
dunque la descrizione interpolatoria e la funzione esatta devono necessariamente coincidere in corrispondenza dei nodi.
Per quanto riguarda una qualche considerazione di convergenza, si può dimostrare che vale la seguente stima per l'errore di interpolazione, nel caso in cui $f \in H^s(-1,1)$:
\[
	\Vert f-\varPi_N^{GLL}f \Vert_{H^k(-1,1)} \leq C_s \left(\frac{1}{N}\right)^{s-k} \Vert f^{(s)} \Vert_{L^2(-1,1)}, \qquad s \geq 1,\ k\in\{0,1\}.
\]
Un risultato di questo tipo è necessario a stabilire la buona posizione dell'impalcatura teorica, poiché garantisce che, a patto di utilizzare un numero di nodi di quadratura sempre maggiore, l'approssimazione discreta della funzione in oggetto tende a quella esatta.\\
Proseguendo con la trattazione, possiamo aggiungere che i pesi della formula di quadratura, definiti da
\[
	\alpha_i = \int_{-1}^1 \psi_i(x)\,\D x,
\]
sono forniti dalla relazione
\begin{equation}
	\alpha_i = \dfrac{2}{N(N+1)} \dfrac{1}{L_N^2(x_i)}.
\end{equation}
In definitiva, la formula di quadratura GLL assume pertanto la forma
\begin{equation}
	\int_{-1}^1 f(x)\,\D x \simeq I_N^{GLL} f := \sum_{k=0}^N \alpha_k f(x_k).
\end{equation}
Tale formula di quadratura ha grado di esattezza pari a $2N-1$, ovvero integra esattamente tutti i polinomi di grado minore o uguale a $2N-1$:
\[
	\int_{-1}^1 f(x)\,\D x = I_N^{GLL}f, \qquad \forall f \in \Q_{2N-1}.
\]
A questo punto è bene procedere a una generalizzazione dei risultati fin qui ottenuti.
Infatti, finora, tutte le analisi sono state svolte prendendo in considerazione solo l'intervallo monodimensionale di riferimento $(-1,1)$.
Invece, poiché intendiamo estendere le formule di quadratura a un generico intervallo $(a,b)$, nodi e pesi andranno definiti come
\[
	t_k = \frac{b-a}{2}x_k + \frac{a+b}{2}, \qquad \delta_k = \frac{b-a}{2} \alpha_k.
\]
In questo modo è immediato definire la formula di quadratura GLL come:
\begin{equation}
	\int_a^b f(x)\,\D x \simeq \sum_{k=0}^N \delta_k f(t_k).
\end{equation}
Altra limitazione che vogliamo assolutamente eliminare è il fatto che, fino a qui, ci siamo occupati esclusivamente di funzioni monodimensionali mentre, in realtà, l'obiettivo finale sarà quello di trattare funzioni a supporto su un dominio bidimensionale.
A questo proposito, l'estensione al caso multidimensionale può essere ottenuta attraverso un prodotto tensoriale della griglia dei nodi GLL fin qui utilizzati.
In pratica, costruiamo una funzione del tipo
\begin{equation}
	\psi_{\vec{k}}(\vec{x}) = \prod_{j=1}^d \psi_{k_j}(x_j),
\end{equation}
dove $d$ indica la dimensione dello spazio nel quale vive la funzione.
In questo modo, la formula di quadratura, nel caso del dominio bidimensionale di riferimento $\hat{\varOmega} = (-1,1)^2$, diventa
\begin{equation}
	\int_{\hat{\varOmega}} f(x_i,x_j)\,\D\varOmega \simeq \sum_{i=0}^N \left(\sum_{j=0}^N f(x_i,x_j)\alpha_j \right)\alpha_i,
\end{equation}
dove i punti $\vec{x}_{ij}=(x_i,x_j)$, con $i,j=0,\dots,N$, rappresentano i nodi GLL su $\hat{\varOmega}$, ottenuti tramite prodotto cartesiano dei nodi GLL sull'intervallo monodimensionale.
Infine, i pesi bidimensionali corrispondono al prodotto dei relativi pesi monodimensionali.
Pertanto avremo
\[
	\alpha_{ij}=\alpha_i\alpha_j.
\]
Gli strumenti fin qui definiti sono sufficienti per l'inquadramento della versione con integrazione numerica del metodo spettrale.
Tuttavia, nel caso in cui si voglia estendere anche il metodo agli elementi spettrali, è necessario generalizzare le formule di integrazione numerica ai singoli elementi $\varOmega_k$, supponendo che essi siano domini quadrangolari di $\R^2$.
Infatti, in questo caso, esiste una trasformazione invertibile $\boldsymbol{\varphi}_k: \hat{\varOmega} \rightarrow \varOmega_k$.\\
I nodi GLL sul generico elemento $\varOmega_k$ vengono generati per mezzo della mappa
\[
	\vec{x}_{ij}^{(k)} = \boldsymbol{\varphi}_k(\vec{x}_{ij}), \qquad i,j=0,\dots,N,
\]
mentre i pesi corrispondenti sono definiti da
\[
	\alpha_{ij}^{(k)} = \alpha_{ij} \vert \det J_k \vert = \alpha_{ij} \dfrac{\vert\varOmega_k\vert}{4}, \qquad i,j=0,\dots,N,
\]
dove con $J_k$ indichiamo lo Jacobiano della trasformazione $\boldsymbol{\varphi}_k$ e con $\vert \varOmega_k \vert$ la misura di $\varOmega_k$.\\
Grazie a queste ultime definizioni, siamo finalmente in grado di fornire la formula GLL su $\varOmega_k$:
\begin{equation}\label{eq:gll}
	\int_{\varOmega_k} f(\vec{x})\,\D\vec{x} \simeq I_{N,k}^{GLL}(f) := \sum_{i,j=0}^N \alpha_{ij}^{(k)}f\left(\vec{x}_{ij}^{(k)}\right).
\end{equation}


\section{I metodi di Galerkin generalizzati}

Una volta introdotti tutti gli elementi necessari a comprendere e implementare l'integrazione numerica, non resta che affrontare l'analisi dei metodi spettrali veri e propri.
Di fatto, le discretizzazioni numeriche di cui andremo a discutere discendono dall'approccio di Galerkin, esattamente come il metodo degli elementi finiti.
Tuttavia, come dicevamo nelle sezioni precedenti, la differenza tra i due approcci riguarda la diversa natura delle funzioni di base che, nel caso spettrale, consistono in polinomi globali, eventualmente anche di alto grado.
Per questo motivo, nonostante la natura completamente diversa dell'approccio utilizzato, parleremo ancora di metodi di Galerkin che però, in questo caso, definiremo generalizzati.

\subsection{Il metodo G-NI}

Il metodo G-NI (\emph{Galerkin with Numerical Integration}) si ottiene approssimando gli integrali all'interno della formulazione SM del problema (\ref{eq:SM}) con le formule di quadratura GLL.
In pratica si sostituisce il prodotto scalare continuo (\ref{eq:prodscal}) con il \emph{prodotto scalare discreto} GLL:
\begin{equation}
	(f,g)_N = \sum_{i,j=0}^N \alpha_{ij}f(\mathbf{x}_{ij})g(\mathbf{x}_{ij}).
\end{equation}
Se ritorniamo al problema (\ref{eq:ellittico}), con $\varOmega = (-1,1)^2$, la sua approssimazione G-NI è data da:
\begin{equation}\label{eq:GNI}
	\mathrm{Trovare}\ u^*_N \in V_N \quad : \quad a_N(u_N^*,v_N)=F_N(v_N), \qquad \forall v_N \in V_N,
\end{equation}
dove $V_N = \{ v \in \Q_N :v \vert_{\varGamma_D} =0\}$, mentre la forma bilineare e il funzionale sono definiti da:
\begin{eqnarray*}
	a_N(u,v)	&=	&(\mu \nabla u,\nabla v)_N + (\sigma u,v)_N \\[2mm]
	F_N(v_N)	&=	&(f,v_N)_N + (h,v_N)_N.
\end{eqnarray*}
Bisogna notare che sia la forma bilineare sia il funzionale non sono più quelli del problema di partenza.
Per questo motivo il metodo G-NI non può essere considerato come un metodo di Galerkin puro bensì, come accennavamo poco prima, come un metodo di Galekin \emph{generalizzato}.\\
Ovviamente i risultati applicabili ai metodi di Galerkin come, ad esempio, il lemma di Ceà, non sono più validi.
Per l'analisi di convergenza del metodo G-NI bisognerà quindi fare riferimento a un risultato diverso.
\begin{lemma}[Strang]
	Si consideri il problema variazionale
	\begin{equation}\label{cont}
		\mathrm{Trovare}\ u \in V \quad : \quad a(u,v)=F(v), \qquad \forall v \in V,
	\end{equation}
	in cui $V$ sia uno spazio di Hilbert con norma $\Vert \cdot \Vert_V,\ F \in V'$ un funzionale lineare e continuo su $V$ e $a(\cdot\,,\cdot):V \times V \rightarrow \mathbb{R}$ una forma bilineare, continua e coerciva su $V$ stesso.
	Sia data inoltre un'approssimazione del problema (\ref{cont}) formulabile attraverso il seguente problema di Galerkin generalizzato:
	\begin{equation}\label{disc}
		\mathrm{Trovare}\ u_N \in V_N \quad : \quad a_N(u_N,v_N)=F_N(v_N), \qquad \forall v_N \in V_,
	\end{equation}
	essendo $\{V_N,N>0\}$ una famiglia di sottospazi di dimensione finita dello spazio infinito dimensionale $V$.\\
	Supponiamo ora che la forma bilineare discreta $a_N(\cdot\,,\cdot)$ sia continua su \mbox{$V_N \times V_N$}, e sia \emph{uniformemente coerciva} su $V_N$, ossia ipotizziamo che \mbox{$\exists\alpha^*>0$}, indipendente da $N$, tale che
	\[
		a_N(v_N,v_N) \geq \alpha^* \Vert v_N \Vert_V^2, \qquad \forall v_N \in V_N.
	\]
	Sia inoltre $F_N$ un funzionale lineare e continuo su $V_N$. Allora:
	\begin{enumerate}
		\item esiste una ed una sola soluzione $u_N$ del problema (\ref{disc});
		\item tale soluzione dipende con continuità dai dati, ossia
		\[
			\Vert u_N \Vert_V \leq \dfrac{1}{\alpha^*} \sup_{v_N \in V_N \backslash \{0\}}\dfrac{F_N(v_N)}{\Vert v_N \Vert_V};
		\]
		\item infine vale la seguente stima a priori dell'errore:
		\begin{eqnarray*}
			\Vert u-u_N \Vert_V	&\leq		&\inf_{w_N \in V_N} \left\{ \left( 1+\dfrac{M}{\alpha^*} \right) \Vert u-w_N \Vert_V\right.\\[2mm]
							&+		&\left.\dfrac{1}{\alpha^*} \sup_{v_N \in V_N \backslash \{0\}} \dfrac{\vert a(w_N,v_N-a_N(w_N,v_N) \vert}{\Vert v_N \Vert_V} \right\}\\[2mm]
							&+		& \dfrac{1}{\alpha^*} \sup_{v_N \in V_N \backslash \{0\}} \dfrac{\vert F(v_N)-F_N(v_N) \vert}{\Vert v_N \Vert_V},
		\end{eqnarray*}
		dove $M$ è la costante di continuità della forma bilineare $a(\cdot\,,\cdot)$.
	\end{enumerate}
\end{lemma}
Per quanto riguarda i dettagli relativi alla dimostrazione, si rimanda a~\cite{QUARTERONI}.\\
A questo punto, invece, non resta che ricavare la formulazione matriciale dell'approssimazione G-NI.
Indichiamo con $\psi_i$, per $i=0,\dots,N$, i polinomi caratteristici associati ai nodi GLL non di Dirichlet.
Tali polinomi costituiscono una base lagrangiana per lo spazio $V_N$, pertanto è possibile scrivere la soluzione come
\[
	u^*_N(x) = \sum_{j=0}^N u^*_N(x_j)\psi_j(x).
\]
Inoltre, possiamo scegliere nella (\ref{eq:GNI}) $v_N = \psi_i$, per $\ i=0,\dots,N$.
In questo modo avremo
\[
	a_N(u^*_N,\psi_i) = F_N(\psi_i), \qquad i=0,\dots,N,
\]
ossia
\[
	\sum_{j=0}^N u^*_N(x_j)a_N(\psi_j,\psi_i) = F_N(\psi_i), \qquad i=0,\dots,N.
\]
Grazie a questi semplici passaggi, la forma matriciale risulta quindi essere:
\begin{equation}\label{matr:gni}
	\A\vec{u}^*_N = \vec{f},
\end{equation}
dove $\vec{u}^*_N = (u^*_N(x_i))$ è il vettore dei coefficienti incogniti, $\A=(a_{ij})$, con $a_{ij}=a_N(\psi_j,\psi_i)$, è la matrice di stiffness associata al problema differenziale.
Invece, per quanto riguarda il vettore $\vec{f}=(f_i)$, definendo $f_i = F_N(\psi_i)$ si mostra banalmente che non è nient'altro se non la discretizzazione del funzionale lineare.\\
In questo modo siamo riusciti a mostrare come la discretizzazione del problema variazionale di partenza abbia condotto alla costruzione di un sistema lineare che però, a differenza di un metodo agli elementi finiti, è caratterizzato da una matrice la cui inversione risulta nettamente più onerosa.

			
\subsection{Il metodo SEM-NI}

Passiamo ora a descrivere la formulazione SEM-NI, ossia il metodo agli elementi spettrali con integrazione gaussiana.
In questo caso il problema variazionale è dato da:
\begin{equation}\label{eq:semni}
	\mathrm{Trovare}\ u_N \in V_N^C\quad : \quad a_{C,N}(u_N,v_N)=F_{C,N}(v_N), \qquad \forall v_N \in V_N^C,
\end{equation}
dove lo spazio $V_N^C$ è ora definito come
\[
	V_N^C = \left\{ v_N \in C^0(\overline{\varOmega}) : v_N\Big\vert_{\varOmega_k} \in \Q_N,\ \forall k=1,\dots,M,\ \mathrm{e}\ v_N\Big\vert_{\varGamma_D} =0\right\},
\]
con $M\in\N$ numero di elementi in cui il dominio $\varOmega$ è stato suddiviso.
Inoltre, definiamo le forme variazionali che compaiono nella formulazione SEM-NI in maniera analoga a quanto fatto in precedenza:
\begin{eqnarray*}
	a_{C,N}(u_N,v_N) &=& \sum_{k} a_{\varOmega_k,N}(u_N,v_N),\\[2mm]
	F_{C,N}(v_N) &=& \sum_k F_{\varOmega_k,N}(v_N).
\end{eqnarray*}	
I termini $a_{\varOmega_k,N}(u_N,v_N)$ e $F_{\varOmega_k,N}(v_N)$ sono ottenuti da
$a_{\varOmega_k}(u_N,v_N)$ e $F_{\varOmega_k}(v_N)$ approssimando ogni integrale con la corrispondente formula di integrazione numerica GLL.\\
A questo punto il nostro obiettivo è trovare una base globale per lo spazio $V_N^C$.
Cominciamo a considerare il caso monodimensionale, con condizioni al bordo di tipo Dirichlet omogenee.
Pertanto scegliamo come dominio in cui ambientare il problema $\hat{\varOmega} = (-1,1)$, partizionandolo in $M$ intervalli disgiunti $\varOmega_m=(\overline{x}_{m-1},\overline{x}_m)$, con $m=1,\dots,M$.
Quando, in un secondo momento, vorremo estendere l'analisi a un contesto multidimensionale, sarà sufficiente, come per il metodo G-NI, aggirare l'ostacolo attraverso un prodotto tensoriale delle funzioni di base.\\
Procediamo dunque con l'analisi del caso unidimensionale.
Innanzitutto si trova una base locale per ogni singolo elemento $\varOmega_m$.
Indicando con $\varphi_m$ la trasformazione affine che mappa l'elemento di riferimento $\hat{\varOmega}$ in $\varOmega_m$, per ogni $m=1,\dots,M$, introduciamo su ogni $\varOmega_m$ l'insieme $\{\psi_i^{(m)}\}_{i=0}^N$ delle funzioni di base tali che
\[
	\psi_i^{(m)}(x) = \psi_i\left(\varphi_m^{-1}(x)\right), \qquad \forall x \in \varOmega_m,
\]
dove $\psi_i$ è il polinomio caratteristico associato al nodo GLL $x_i$ in $\hat{\varOmega}$.\\
La soluzione discreta $u_N$ del SEM, su ciascun sotto-elemento, può essere definita come
\begin{equation}
	u_N(x) = \sum_{i=0}^N u_i^{(m)} \psi_i^{(m)}(x),
\end{equation}
con $u_i^{(m)}=u_N(x_i^{(m)})$ i coefficienti dello sviluppo in serie, dove $x_i^{(m)}$ non è altro che l'immagine attraverso la mappa $\varphi_m$ \mbox{dell'$i$-esimo} nodo GLL di $\hat{\varOmega}$.\\
Le funzioni di base associate ai nodi interni di $\varOmega_m$ assumono quindi la forma
\[
	\tilde{\psi}_i^{(m)}(x) =
	\left\{
		\begin{array}{ll}
			\psi_i^{(m)}(x)		&\qquad \textrm{per\ } x \in \varOmega_m\\[2mm]
			0				&\qquad \textrm{altrimenti}.
		\end{array}
	\right.
\]
Mentre, per quanto riguarda i nodi estremi $\overline{x}_m$ dei sottodomini $\varOmega_m$, con $m=1,\dots,M-1$, si definiscono invece le funzioni di base
\[
	\psi_m^*(x) =
	\left\{
		\begin{array}{ll}
			\psi_N^{(m)}(x)		&\qquad\textrm{per\ } x \in \varOmega_m\\[2mm]
			\psi_0^{(m+1)}(x)	&\qquad\textrm{per\ } x \in \varOmega_{m+1}\\[3mm]
			0				&\qquad\textrm{altrimenti}.
		\end{array}
	\right.
\]
È bene sottolineare che le funzioni di base, in corrispondenza dei nodi estremi $\overline{x}_0$ e $\overline{x}_M$, non sono definite, in quanto in tali punti sono state poste condizioni Dirichlet.
Invece, ad esempio, nel caso in cui avessimo a che fare con condizioni di tipo Neumann, esisterebbero sia $\psi_0^* = \psi_0^{(1)}$ sia $\psi_M^* = \psi_N^{(M)}$.\\
A questo punto, avendo definito una base globale per lo spazio $V_N^C$, possiamo esprimere ogni funzione $u_N \in V_N^C$ semplicemente come
\[
	u_N(x) = \sum_{m=1}^{M-1} u_m^{\varGamma} \psi_m^*(x) + \sum_{m=1}^M \sum_{i=1}^{N-1} u_i^{(m)} \tilde{\psi}_i^{(m)}(x),
\]
dove $u_m^{\varGamma} = u_N(\overline{x}_m)$ sono i coefficienti relativi ai nodi di bordo di ogni singolo sottodominio.


\section{Problemi di evoluzione}

Dopo aver introdotto tutto l'apparato teorico necessario alla risoluzione di problemi ellittici attraverso l'approccio spettrale, finalmente veniamo a occuparci della parte sulla quale è incentrato il vero obiettivo di questa analisi: la risoluzione di problemi tempo dipendenti.
Di fatto, la difficoltà nell'affrontare la risoluzione di problemi evolutivi risiede unicamente negli schemi iterativi che consentono di ottenere la soluzione ad ogni istante temporale.
L'obiettivo del lavoro sarà dunque quello di implementare un solutore per l'avanzamento temporale, accoppiandolo con un risolutore spettrale per il problema ellittico che dovrà essere risolto a ogni passo.
In particolare, in questo lavoro abbiamo focalizzato la nostra attenzione su quella classe di risolutori che va sotto il nome di $\vartheta$-metodo.
Si tratta di un approccio alle differenze finite classico ad un passo che però, grazie all'introduzione di un parametro di rilassamento, permette di ottenere una vasta gamma di metodi, sia impliciti sia espliciti.
Prima di procedere con la descrizione dell'algoritmo e della sua analisi, iniziamo con l'inquadrare dal punto di vista analitico i problemi evolutivi che intendiamo risolvere in questo lavoro: i problemi parabolici.

\subsection{Problemi parabolici}

Sia $\varOmega \subset \R^2$ un dominio (aperto connesso) limitato, $T\in(0,+\infty)$ un numero strettamente positivo, e consideriamo il cilindro spazio-temporale $Q_T = \varOmega \times (0,T)$.\\
Una volta assegnata la forzante $f$ in $Q_T$, vogliamo determinare una soluzione $u$ dell'equazione parabolica
\begin{equation}\label{eq:parabolic}
	\Deriv{}{u}{t} + Lu = f, \qquad \textrm{in\ } Q_T
\end{equation}
che soddisfi la condizione iniziale di Cauchy
\[
	u(\vec{x},0) = u_0(\vec{x})
\]
unitamente alle condizioni al contorno di tipo Dirichlet, Neumann, Robin o miste sulla frontiera laterale di $Q_T$
\[
	S_T=\partial \varOmega \times [0,T],
\]
come, ad esempio,
\begin{equation}
	\left\{
	\begin{array}{ll}
		u(\vec{x},t) = \varphi(\vec{x},t)				&\qquad\textrm{per\ }\vec{x}\in\varGamma_D,\quad t>0\\[3mm]
		\nabla u(\vec{x},t)\cdot\vec{n} = \psi(\vec{x},t)	&\qquad\textrm{per\ }\vec{x}\in\varGamma_N,\quad t>0
	\end{array}
	\right.
\end{equation}
dove $\varphi$ e $\psi$ sono funzioni assegnate, e $\varGamma_D \cup \varGamma_N = \partial \varOmega$, con $\varGamma_D \cap \varGamma_N = \emptyset$.
Naturalmente, come consuetudine, rappresentiamo la porzione di bordo sulla quale intendiamo imporre condizioni di tipo Dirichlet con $\varGamma_D$, mentre quella sulla quale intendiamo imporre condizioni di tipo Neumann verrà indicata con $\varGamma_N$.\\
Da ora in poi consideriamo un problema test, che verrà inquadrato dal punto di vista teorico, in modo da mostrare le proprietà caratterizzanti la classe dei problemi evolutivi.
Supporremo che l'operatore ellittico sia del secondo ordine, in particolare che $L=-\varDelta$, così da analizzare la classe di problemi che riguarda la diffusione del calore in un mezzo.

A tal proposito, cerchiamo un'ambientazione funzionale appropriata.\\
Per quanto riguarda la forzante $f$, una condizione naturale da porre è che sia $f \in L^2(Q_T)$.
Questo equivale a garantire che
\[
	f \in L^2(0,T;L^2(\varOmega)).
\]
Coerentemente, la condizione di Dirichlet si traduce, pensando alla formulazione debole del problema, nella richiesta che $u$ appartenga allo spazio $V=H^1_{\varGamma_D}(\varOmega)$ per quasi ogni $t \in [0,T]$.
Inoltre, si richiede che la funzione $t \mapsto \Vert u(t) \Vert_V$ appartenga a $L^2(0,T)$.
In questo modo intendiamo imporre che
\[
	u \in L^2(0,T;V).
\]
Infine, la richiesta non troppo restrittiva sul dato iniziale $u_0(\vec{x})$ è che
\[
	u_0 \in L^2(\varOmega).
\]
Consideriamo ora il termine $\partial u / \partial t$.
Dato che $u(t) \in V$, si ha che $\varDelta u(t) \in V'$, cioè $H^{-1}(\varOmega)$ per quasi ogni $t>0$.
Per questo motivo, la condizione naturale da imporre è che $\partial u / \partial t \in V'$ per quasi ogni $t>0$.
Analogamente, sempre per questioni riguardanti la concordanza tra spazi funzionali, si richiede che la funzione $t \mapsto \Vert \partial u(t) / \partial t \Vert_{V'}$ appartenga a $L^2(0,T)$.
Quindi, in definitiva, possiamo unire tutte queste considerazioni imponendo che
\[
	\Deriv{}{u}{t} \in L^2(0,T;V').
\]
Grazie a tutte queste premesse, e sfruttando alcuni risultati di analisi funzionale, per i cui dettagli rimandiamo a~\cite{SALSA}, è possibile dedurre che
\[
	u \in C^0([0,T];L^2(\varOmega)).
\]
Fino ad ora abbiamo, di fatto, introdotto degli ambienti funzionali spazio-temporali, che sono necessari a descrivere correttamente i problemi parabolici.
In realtà, ognuno di questi spazi può essere facilmente descritto grazie alle definizioni:
\[
	\begin{array}{c}
		L^2([0,T];W) = \left\{v:[0,T] \rightarrow W: t \rightarrow v(t)\ \mathrm{misurabile\ e}\ \displaystyle{\int_0^T \Vert v(t) \Vert_W^2} < \infty \right\}\\[6mm]
		C^0([0,T];W) = \left\{v:[0,T] \rightarrow W: \Vert v(t) \Vert_W \in C^0([0,T]), \forall t \geq 0\right\},
	\end{array}
\]
dove $W$ è uno spazio di Hilbert con norma $\Vert \cdot \Vert_W$ che, a seconda dei casi, può trattarsi dello spazio $V$ appena definito, oppure del suo duale.\\
Dopo aver introdotto tutte gli ambienti funzionali corretti, in modo da garantire l'esistenza della soluzione, passiamo a descrivere la formulazione debole del problema:
\begin{equation}
	\int_{\varOmega} \Deriv{}{u(t)}{t}v\,\D\varOmega + a(u(t),v)=F(v), \qquad \forall v \in V,
\end{equation}
con $u(0)=u_0$, dove $a(\cdot\,,\cdot)$ e $F(\cdot)$ sono la forma bilineare e il funzionale associati, rispettivamente, all'operatore ellittico e al termine noto.
Nel caso di problemi parabolici, la condizione di coercività della forma bilineare può essere rilassata, garantendo comunque la buona posizione del problema.
La condizione che richiediamo in questo caso è, pertanto, la \emph{coercività debole}, ovverosia che
\[
	\exists \lambda \geq 0,\ \exists \alpha >0 \quad : \quad a(v,v)+\lambda \Vert v \Vert_{L^2(\varOmega)}^2 \geq \alpha \Vert v \Vert_V^2, \qquad \forall v \in V.
\]
Ricordiamo che, una volta garantita questa proprietà, è sempre possibile ricondursi a una forma coerciva in senso stretto, semplicemente attraverso un cambio di variabili.
D'altra parte, come nel caso dei problemi ellittici, ricordiamo che è sempre necessario dimostrare la linearità e continuità del funzionale, oltre alla bilinearità e continuità della forma $a(\cdot\,,\cdot)$.
Sotto queste ipotesi, esistenza, unicità e dipendenza con continuità dai dati possono essere dimostrate attraverso l'uso del metodo di Faedo-Galerkin.
Anche in questo caso, se volessimo trattare l'argomento in maniera approfondita, svieremmo troppo dal discorso principale.
Per questo motivo rimandiamo ogni dettaglio riguardante questi aspetti teorici a~\cite{SALSA}.

Una volta definita la formulazione debole, ambientando ogni termine nel corretto spazio funzionale, l'approssimazione G-NI del problema (\ref{eq:parabolic}) si ottiene tramite la semi-discretizzazione spaziale della formulazione variazionale.
In pratica, \mbox{$\forall t>0$}, possiamo definirla come
\begin{equation}
	\begin{array}{c}
		\mathrm{Trovare}\ u^*_N(t) \in V_N:\\[3mm]
		\displaystyle{\int_{\varOmega}} \Deriv{}{u^*_N(t)}{t}v_N\,\D\varOmega + a_N\left(u^*_N(t),v_N\right) = F_N(v_N), \qquad \forall v_N \in V_N,
	\end{array}
\end{equation}
con $u^*_N(0) = u_{0N}$, dove $u_{0N}$ è la proiezione di $u_0$ nello spazio $V_N$, approssimazione finito dimensionale dello spazio funzionale $V$.
Per studiare le proprietà di stabilità di questo problema, come è solito fare, scegliamo come funzione test proprio l'approssimazione numerica della soluzione, e otteniamo
\[
	\int_{\varOmega} \Deriv{}{u^*_N(t)}{t}u^*_N(t)\,\D\varOmega + a_N\left(u^*_N(t),u^*_N(t)\right) = F_N(u^*_N(t)).
\]
A questo punto, sfruttando la coercività debole e la definizione di derivata di un prodotto, possiamo ricavare
\[
	\dfrac{1}{2}\Deriv{}{}{t}\Vert u^*_N(t) \Vert_{L^2(\varOmega)}^2 + \alpha \Vert u^*_N(t) \Vert_V^2 \leq F_N(u^*_N(t)) + \lambda\Vert u^*_N(t) \Vert_{L^2(\varOmega)}^2
\]
dove, come già osservato in precedenza, $\alpha$ è la costante di coercività per la forma bilineare.
Ora, applicando la disuguaglianza di Young al secondo termine, con $\varepsilon=\alpha/2$, possiamo trovare una costante $C$, che dipenda solo da $\alpha$ ma non da $N$, in modo che $\forall t>0$
\[
	\Vert u^*_N(t) \Vert_{L^2(\varOmega)}^2 + \alpha\int_0^T\Vert u^*_N(t) \Vert_V^2\,\D t \leq \Vert u_{0N} \Vert_{L^2(\varOmega)}^2 + C\int_0^T\Vert f(t) \Vert^2\,\D t.
\]
In questo modo siamo riusciti a dimostrare la stabilità dell'approssimazione di Galerkin nel caso spettrale.
Incidentalmente notiamo che non possiamo svolgere un'analisi simile a quella che si potrebbe fare nel caso di elementi finiti poiché, visto che ci stiamo occupando di approssimazioni G-NI, dobbiamo fondare tutte le nostre affermazioni solo sul lemma di Strang.\\
In ogni caso, grazie ai fondamenti teorici, siamo anche in grado di fornire qualche risultato sulla convergenza dell'approssimazione che stiamo analizzando.
Introduciamo la variabile $e(t) = R_Nu(t)-u^*_N(t)$, dove $R_N$ è l'operatore di proiezione che mappa la soluzione esatta nello spazio discreto che stiamo considerando.
La funzione di errore $e(t)$, allora, potrà soddisfare la relazione
\[
	\dfrac{1}{2}\Deriv{}{}{t}\Vert e(t) \Vert_{L^2(\varOmega)}^2 + \alpha \Vert e(t) \Vert_V^2 \leq \left\vert\left(u_t - R_N u_t, e(t)\right) + a_N\left(u-R_Nu, e(t)\right)\right\vert.
\]
A questo punto, grazie al teorema di rappresentazione di Riesz, per $u\in\ V$ possiamo definire la nuova norma
\[
	\Vert u \Vert_{V'} = \sup_{v\in V, v\neq 0} \dfrac{(u,v)}{\Vert v \Vert_V}
\]
che è semplicemente la norma di $u$ nello spazio duale $V'$.
In questo modo, usando la definizione appena enunciata insieme alla continuità dell'operatore ellittico, se ne deduce che
\[
	\begin{array}{c}
		\left\vert\left(u_t - R_N u_t, e(t)\right) + a_N\left(u-R_Nu, e(t)\right)\right\vert \leq\\[3mm]
		\leq C\Vert e(t) \Vert_V \bigg[\Vert u_t - R_N u_t \Vert_{V'} + \Vert u-R_N u \Vert_V \bigg].
	\end{array}
\]
Grazie a quest'ultima relazione è ora immediato ricavare la stima
\[
	\begin{array}{c}
		\Vert e(t) \Vert_{L^2(\varOmega)}^2 + \alpha\displaystyle{\int_0^T\Vert e(t) \Vert_V^2\,\D t}\leq\\[3mm]
		\leq \Vert e(0) \Vert_{L^2(\varOmega)}^2 + C \left[ \displaystyle{\int_0^T \Vert(u_t-R_Nu_t)(t)\Vert_{V'}^2\,\D t} + \displaystyle{\int_0^T \Vert (u-R_Nu)(t)\Vert_{V}^2\,\D t}\right],
	\end{array}
\]
dove $C$ è una costante indipendente da $N$.
Pertanto, possiamo concludere che lo schema proposto è convergente se ciascun addendo nel secondo termine tende a $0$, quando $N\to\infty$, sotto l'ipotesi che $u$, $u_t$ e $u_0$ siano sufficientemente regolari.
In ogni caso, come accade sempre nei problemi parabolici, la regolarità della soluzione segue direttamente da quella dei dati iniziali e al bordo, oltre che dal soddisfacimento di alcune condizioni di compatibilità tra loro, di cui non daremo ulteriori dettagli.
Ciò che interessava, a questo livello della discussione, era unicamente mostrare come, anche nel caso di approssimazione spettrale G-NI, sia sempre possibile ricavare alcuni risultati di stabilità e convergenza.\\
A questo punto, la formulazione algebrica può essere ottenuta procedendo in maniera analoga a quanto fatto per il metodo G-NI nel caso ellittico.
Di fatto, si sceglie $\{\psi_j\}$ come base dello spazio $V_N$ e, per ogni $t>0$, si scrive la soluzione secondo lo sviluppo in serie
\[
	u^*_N(\vec{x},t) = \sum_{j=0}^N u_j(t) \psi_j(\vec{x}).
\]
Ora, a differenza del caso ellittico, i coefficienti $\{u_j(t)\}$, che rappresentano le incognite del problema, sono dipendenti dall'istante temporale $t$.\\
Indicando con $u'_j(t)$ la derivata della funzione $u_j(t)$ rispetto a $t$, si ottiene
\[
	\int_{\varOmega} \sum_{j=0}^N u'_j(t) \psi_j \psi_i\,\D\varOmega + a_N \left(\sum_{j=0}^N u_j(t) \psi_j,\psi_i \right) = F_N(\psi_i), \qquad i=0,\dots,N,
\]
ossia, scambiando in maniera opportuna gli operatori di integrazione con le serie,
\[
	\sum_{j=0}^N u'_j(t) \underbrace{\int_{\varOmega} \psi_j \psi_i\,\D\varOmega}_{m_{ij}} + \sum_{j=0}^N u_j(t) \underbrace{a_N(\psi_j,\psi_i)}_{a_{ij}} = \underbrace{F_N(\psi_i)}_{f_i(t)}, \qquad i=0,\dots,N.
\]
Dopo aver definito le varie forme lineari e bilineari nella precedente relazione, è ora possibile ottenere il vettore delle incognite $\vec{u} = [u_1(t),\dots,u_N(t)]^T$, la matrice di massa $\M=[m_{ij}]$, la matrice di stiffness $\A=[a_{ij}]$ e il vettore dei termini noti $\vec{f}=[f_1(t),\dots,f_N(t)]^T$.
In questo modo, grazie all'introduzione del formalismo algebrico, il sistema può essere riscritto in maniera compatta:
\begin{equation}
	\M \vec{u}'(t) + \A\vec{u}(t) = \vec{f}(t).
\end{equation}
Abbiamo pertanto ottenuto il sistema lineare, tempo dipendente, che dovremo risolvere per trovare la soluzione numerica del problema parabolico.
Di fatto, a questo punto, abbiamo a che fare con un sistema di equazioni differenziali ordinarie, poiché in effetti l'incognita è vettoriale.
Con l'obiettivo di risolvere sistemi di questo tipo, passiamo a mostrare i metodi numerici che intendiamo applicare in questo lavoro.

\subsection{Il $\vartheta$-metodo}

Come dicevamo nell'introduzione a questa sezione, la discretizzazione temporale può essere effettuata con il cosiddetto \emph{$\vartheta$-metodo}, che consiste nel discretizzare la derivata temporale con un rapporto incrementale, sostituendo gli altri termini con una combinazione lineare convessa, dipendente dal parametro reale $\vartheta\in[0,1]$, del valore al tempo $t^k$ e di quello al tempo $t^{k+1}$:
\begin{equation}\label{eq:thetametodo}
	\M \dfrac{\vec{u}^{k+1}-\vec{u}^k}{\Delta t}+\A[\vartheta\vec{u}^{k+1}+(1-\vartheta)\vec{u}^k]=\vartheta \vec{f}^{k+1}+(1-\vartheta)\vec{f}^k.
\end{equation}
Il passo di discretizzazione temporale (supposto costante) è indicato con \mbox{$\Delta t=t^{k+1}-t^k$}, per $k=0,1,\dots$, mentre il sopraindice $k$ indica che la quantità in considerazione è riferita al tempo $t^k$.
Per particolari valori di $\vartheta$ si ottengono metodi di discretizzazione temporale ben noti:
\begin{itemize}
	\item per $\vartheta=0$ si ottiene il metodo di \emph{Eulero esplicito}:
	\[
		\M \dfrac{\vec{u}^{k+1}-\vec{u}^k}{\Delta t}+\A\vec{u}^k = \vec{f}^k,
	\]
	\item per $\vartheta=1$ si ottiene il metodo di \emph{Eulero implicito}:
	\[
		\M \dfrac{\vec{u}^{k+1}-\vec{u}^k}{\Delta t}+\A\vec{u}^{k+1} = \vec{f}^{k+1},
	\]
	\item mentre per $\vartheta=1/2$ si ottiene il metodo di \emph{Crank-Nicolson}:
	\[
		\M \dfrac{\vec{u}^{k+1}-\vec{u}^k}{\Delta t}+\dfrac{1}{2}\A(\vec{u}^{k+1}+\vec{u}^k)=\dfrac{1}{2}(\vec{f}^{k+1}+\vec{f}^k).
	\]
\end{itemize}
Il $\vartheta$-metodo è accurato al primo ordine, ad eccezione del caso $\vartheta=1/2$, la cui accuratezza è del secondo ordine.
Inoltre è importante sottolineare che il metodo risulta essere esplicito per il solo valore $\vartheta=0$ mentre, in tutti gli altri casi, lo schema è sempre implicito.\\
Per quanto riguarda la stabilità del metodo di discretizzazione temporale, si può dimostrare (si veda a tal proposito~\cite{QUARTERONI}) che, se $\vartheta \geq 1/2$, allora il $\vartheta$-metodo è incondizionatamente stabile, ossia è stabile per ogni $\Delta t$.
Invece, se $\vartheta <1/2$, si ha assoluta stabilità solo se $\Delta t \leq C(\vartheta)N^{-4}$, dove $C(\vartheta)$ è una costante positiva, monotona crescente in $\vartheta$.\\

\subsection{L'equazione del calore}

Passiamo ora a descrivere alcuni risultati di stabilità e convergenza per un caso test di equazione del calore che, per semplicità, supponiamo monodimensionale.
Nonostante la trattazione sia decisamente più complicata rispetto al caso in cui il problema è discretizzato attraverso l'impiego di elementi finiti, i risultati possono comunque essere estesi anche al caso multidimensionale.\\
Consideriamo allora l'equazione del calore con condizioni al bordo di Dirichlet omogenee:
\begin{equation}
	\left\{
	\begin{array}{ll}
		u_t - u_{xx} = 0,		&\qquad -1<x<1,\quad t>0\\[2mm]
		u(x,0) = u_0(x),		&\qquad -1<x<1\\[2mm]
		u(-1,t) = u(1,t) = 0,	&\qquad t>0.
	\end{array}
	\right.
\end{equation}
Applichiamo la discretizzazione spaziale G-NI e la discretizzazione temporale attraverso il $\vartheta$-metodo.\\
Siano $\Delta t>0$ il passo di discretizzazione temporale e $t^k = k\Delta t$.
In questo modo la discretizzazione completa dell'equazione del calore diventa:
\begin{equation}\label{eq:heatspect}
	\begin{array}{c}
		\forall k \geq 0, \mathrm{trovare}\ u^{N,k}:\\[3mm]
		\left(u^{N,k+1}-u^{N,k},v\right)_N + \Delta t\,a_N\left(\vartheta u^{N,k+1}+(1-\vartheta)u^{N,k},v\right) = 0.
	\end{array}
\end{equation}
Inoltre è necessario introdurre $u^{N,0}$, ovverosia l'interpolante di $u_0$ sui nodi GLL, indispensabile per trattare il dato iniziale in maniera coerente rispetto alla formulazione appena definita.
Da ora in poi, per semplicità di notazione, indichiamo con $\Vert v \Vert_0$ la norma $L^2(-1,1)$ di $v$.\\
Per provare la stabilità del metodo, scegliamo $v = \vartheta u^{N,k+1}+(1-\vartheta)u^{N,k}$.
In questo modo otteniamo che
\begin{equation}
	\begin{array}{ll}
		\vartheta \Vert u^{N,k+1} \Vert_N^2	& + (1-2\vartheta) (u^{N,k+1},u^{N,k})_N+\\[4mm]
									& - (1-\vartheta)\Vert u^{N,k} \Vert_N^2+\\[2mm]
									& +\dfrac{\Delta t}{4} \left\Vert \vartheta u_x^{N,k+1} + (1-\vartheta)u_x^{N,k} \right\Vert_0^2 \leq 0.
	\end{array}
\end{equation}
Ora, dato che $1-2\vartheta \leq 0$, grazie alla disuguaglianza di Cauchy-Schwartz possiamo dedurre che
\begin{equation}
	(1-2\vartheta) \left(u^{N,k+1},u^{N,k}\right)_N \geq \left(\dfrac{1}{2}-\vartheta\right)\left(\Vert u^{N,k+1} \Vert_N^2 + \Vert u^{N,k} \Vert_N^2 \right)
\end{equation}
e, di conseguenza, che
\begin{equation}
	\Vert u^{N,k+1} \Vert_N^2 + \dfrac{\Delta t}{2} \left\Vert \vartheta u_x^{N,k+1} +(1-\vartheta)u_x^{N,k} \right\Vert_0^2 \leq \Vert u^{N,k} \Vert_N^2.
\end{equation}
In questo modo, quindi, per ogni $k \geq 0$ ricaviamo immediatamente che
\begin{equation}
	\Vert u^{N,k} \Vert_N^2 + \dfrac{\Delta t}{2} \sum_{j=0}^{k-1} \left\Vert \vartheta u_x^{N,j+1}+ (1-\vartheta) u_x^{N,j} \right\Vert_0^2 \leq \Vert u_0 \Vert_N^2.
\end{equation}
Grazie a quest'ultima relazione possiamo notare che lo schema è incondizionatamente stabile se $\vartheta \in [ \frac{1}{2},1 ]$.
Inoltre, riusciamo anche a mostrare che lo schema è contrattivo, questo poiché $\Vert u^{N,k+1} \Vert_N \leq \Vert u^{N,k} \Vert_N$, per ogni $k \geq 0$.

Passiamo ora a verificare la convergenza del metodo.\\
A tal proposito definiamo la proiezione $\tilde{u}(t) = \varPi_N u(t) \in \mathbb{P}_N^0$.
In questo modo, per ogni $v \in \mathbb{P}_N^0$ si ha che
\begin{equation}
	\left(\vartheta \tilde{u}_t^{k+1} + (1-\vartheta)\tilde{u}_t^k,v\right) + a\left(\vartheta \tilde{u}^{k+1} + (1-\vartheta)\tilde{u}^k,v\right) = (\delta^k,v),
\end{equation}
dove $\delta^k = \vartheta (\tilde{u}-u)_t^{k+1}+(1-\vartheta)(\tilde{u}-u)_t^k$.
Definiamo ora l'errore $e^k = u^{N,k}-\tilde{u}^k$.
Sfruttando la relazione (\ref{eq:heatspect}) si ottiene subito che
\begin{equation}
	\begin{array}{ll}
		\dfrac{1}{\Delta t}\left(u^{N,k+1}	-u^{N,k},v\right)_N	& -\left(\vartheta \tilde{u}_t^{k+1}+(1-\vartheta)\tilde{u}_t^k,v\right)_N+\\[4mm]
												& + a_N\left(\vartheta e^{k+1}+(1-\vartheta)e^k,v\right)=\\[4mm]
												& = -(\delta^k,v) - E\left(\vartheta \tilde{u}_t^{k+1} + (1-\vartheta)\tilde{u}_t^k,v\right),
	\end{array}
\end{equation}
dove la forma bilineare $E$ è definita come $E(\varphi,\psi) = (\varphi,\psi)-(\varphi,\psi)_N$.\\
Sia ora $z = z(t)$ una funzione differenziabile con continuità sull'intervallo $(0,+\infty)$.
In questo modo possiamo definire
\[
	\epsilon^k(z) = \dfrac{1}{\Delta t}(z^{k+1}-z^k)-\left(\vartheta z_t^{k+1}+(1-\vartheta)z_t^k\right).
\]
Se $z \in C^2(0,+\infty)$, applicando la formula di Taylor otteniamo che
\[
	\epsilon^k(z) = \dfrac{1}{\Delta t} \int_{t^k}^{t^{k+1}} \left(s-(1-\vartheta)t^{k+1}- \vartheta t^k\right) z_{tt}(s)\,\D s,
\]
da cui si ricava immediatamente che
\begin{equation}\label{taylor:theta}
	\vert \epsilon^k(z) \vert \leq \max (\vartheta,1-\vartheta) \int_{t^k}^{t^{k+1}} \vert z_{tt}(s) \vert\,\D s \leq \int_{t^k}^{t^{k+1}} \vert z_{tt}(s)\vert\,\D s.
\end{equation}
A questo punto, se $\vartheta =1/2$ e $z \in C^3(0,+\infty)$, allora è possibile ottenere una stima migliore attraverso una formula di Taylor di ordine più elevato, cioè
\[
	\epsilon^k(z) = \dfrac{1}{2\Delta t} \int_{t^k}^{t^{k+1}} \left(t^k-s\right)\left(t^{k+1}-s\right) z_{ttt}(s)\,\D s,
\]
da cui possiamo dedurre che
\begin{equation}\label{taylor:cn}
	\vert \epsilon^k(z) \vert \leq \dfrac{\Delta t}{8} \int_{t^k}^{t^{k+1}} \vert z_{ttt}(s)\vert\,\D s.
\end{equation}
Ora, scegliendo $v = \vartheta e^{k+1}+(1-\vartheta)e^k$ e sfruttando alcuni risultati di analisi funzionale come la disuguaglianza di Cauchy-Schwartz, la disuguaglianza di Poincaré, la disuguaglianza triangolare e il teorema di equivalenza fra norma discreta e norma continua (per i cui dettagli rimandiamo a~\cite{SPETTRALI}), dopo diversi passaggi si ottengono le stime di convergenza nelle norme $L^2$ e $H^1$.\\
In particolare, nel caso in cui $\vartheta =1/2$, per la norma $L^2$ ricaviamo che
\begin{equation}
	\begin{array}{ll}
		\left\Vert u(t^k)-u^{N,k} \right\Vert_{L^2}\leq	& C_1 N^{-r} \left(\vert u_0 \vert_{H^r}^2 + \vert u(t^k) \vert_{H^r}^2 + \Delta t\, \displaystyle{\sum_{j=0}^k} \vert u_t^j \vert_{H^r}^2 \right)^{1/2} \\[6mm]
											& + C_2 \Delta t^2 \left(\displaystyle{\int_0^{t^k} \Vert \tilde{u}_{ttt}(s)} \Vert_{L^2}^2\,\D s \right)^{1/2}
	\end{array}
\end{equation}
mentre, per quando riguarda la norma $H^1$, otteniamo
\begin{equation}
	\begin{array}{c}
		\left( \Delta t\, \displaystyle{\sum_{j=0}^{k-1}} \left\Vert \vartheta \left(u(t^j)-u^{N,j}\right) + (1-\vartheta)\left(u(t^{j+1})-u^{N,j+1}\right) \right\Vert_{H^1}^2 \right)^{1/2} \\[6mm]
		\leq C_3 N^{1-r} \left( \vert u_0 \vert_{H^{r-1}}^2 + \vert u(t^k) \vert_{H^r}^2 + \Delta t\,\displaystyle{\sum_{j=0}^k} \left\vert u_t^j \right\vert_{H^{r-1}}^2 \right)^{1/2} \\[6mm]
		+ C_2 \Delta t^2 \left(\displaystyle{\int_0^{t^k} \Vert \tilde{u}_{ttt}(s) \Vert_{L^2}^2\,\D s} \right)^{1/2}.
	\end{array}
\end{equation}
Infine, ricordando che $\tilde{u}_t = \varPi_N u_t$, si può sostituire la derivata temporale di $\tilde{u}$ con quella di $u$.\\
Concludendo il discorso, ricordiamo che queste stime di convergenza valgono solo nel caso in cui $\vartheta=1/2$.
Per quanto concerne lo scenario in cui $\vartheta \in (\frac{1}{2},1]$, utilizzando la (\ref{taylor:theta}) invece della (\ref{taylor:cn}), si possono ottenere stime analoghe, con l'unica differenza che figurerà $\Delta t$ al posto di $\Delta t^2$ e $u_{tt}$ invece di $u_{ttt}$.
